
### Creates figure 1
### Prevalence of headscarf and burqa laws


# Set up ------------------------------------------------------------------

# set working directory
setwd("")

# load packages
packs <- c("foreign", "ggplot2", "ggthemes", "WDI", "countrycode", "directlabels", "maptools", 
           "ggmap", "sp", "grid", "Hmisc", "cowplot", "rgdal", "grid", "gridExtra")
lapply(packs, require, character.only = TRUE)


# Preparing data for maps -----------------------------------------

# read in world shapefile
countries <- readOGR(dsn = "ne_50m_admin_0_countries.shp",
                     layer="ne_50m_admin_0_countries")

# keep europe and bordering asian countries
countries <- subset(countries, countries$CONTINENT=="Europe" | countries$CONTINENT=="Asia")

# fortify
countries_df<-fortify(countries, region="ADM0_A3_IS")

# load veiling data
veil <- read.csv("nuts0_bans.csv")

# merge veiling and world map
veil_map <- merge(countries_df, veil, by="id", all.x=T)

# order the dataset
veil_map <- veil_map[order(veil_map$order),]


# Building maps -----------------------------------------

## Headscarf ban map
hscrf.map <- ggplot() +
  geom_polygon(data=veil_map, aes(x=long, y=lat, group=group, fill=as.factor(hscrf_ban)), 
               alpha=.4) +
  geom_path(data=veil_map, aes(long,lat, group=group), color="black",
            size=0.3) +    
  coord_cartesian(xlim = c(-21, 42), ylim=c(35, 70)) +
  scale_fill_grey(start = 0, end = 0.9, labels=c("no ban", "local", "national"), 
                  name="", na.translate=FALSE) +
  labs(title = "Headscarf Ban") +
  theme_bw() +
  theme(aspect.ratio=.8,
        axis.title=element_blank(),
        axis.text=element_blank(),
        axis.ticks=element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        plot.title = element_text(size=16, face="bold", hjust=.5),
        legend.text=element_text(size=12),
        panel.background = element_rect(fill = 'white', colour = NA),
        legend.position="bottom")

## Burqa ban map
burqa.map <- ggplot() +
  geom_polygon(data=veil_map, aes(x=long, y=lat, group=group, fill=as.factor(burqa_ban)), 
               alpha=.4) +
  geom_path(data=veil_map, aes(long,lat, group=group), color="black",
            size=0.3) +    
  coord_cartesian(xlim = c(-21, 42), ylim=c(35, 70)) +
  scale_fill_grey(start = 0, end = 0.9, labels=c("no ban", "local", "national"), 
                  name="", na.translate=FALSE) +
  labs(title = "Burqa Ban") +
  theme_bw() +
  theme(aspect.ratio=.8,
        axis.title=element_blank(),
        axis.text=element_blank(),
        axis.ticks=element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        plot.title = element_text(size=16, face="bold", hjust=.5),
        legend.text=element_text(size=12),
        panel.background = element_rect(fill = 'white', colour = NA),
        legend.position="bottom")



# save maps
ggsave("hscrf_map.pdf", hscrf.map)
ggsave("burqa_map.pdf", burqa.map)


